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ABSTRACT 



Context. The poorly-ionized interior of the protoplanetary disk or 'dead zone' is the location where dust coagulation processes may 
be most efficient. However even here, planetesimal formation may be limited by the loss of solid material through radial drift, and by 
collisional fragmentation of the particles. Both depend on the turbulent properties of the gas. 

Aims. Our aim here is to investigate the possibility that solid particles are trapped at local pressure maxima in the dynamically 
evolving disk. We perform the first 3-D global non-ideal magnetohydrodynamical (MHD) calculations of a section of the disk treating 
the turbulence driven by the magneto-rotational instability (MRI). 

Methods. We use the ZeusMP code with a fixed Ohmic resistivity distribution. The domain contains an inner MRI-active region near 
the young star and an outer midplane dead zone, with the transition between the two modeled by a sharp increase in the magnetic 
diffusivity. 

Results. The azimuthal magnetic fields generated in the active zone oscillate over time, changing sign about every 150 years. We 
thus observe the radial structure of the 'butterfly pattern' seen previously in local shearing-box simulations. The mean magnetic field 
diffuses from the active zone into the dead zone, where the Reynolds stress nevertheless dominates, giving a residual a between 10 -4 
and 10~ 3 . The greater total accretion stress in the active zone leads to a net reduction in the surface density, so that after 800 years an 
approximate steady state is reached in which a local radial maximum in the midplane pressure lies near the transition radius. We also 
observe the formation of density ridges within the active zone. 

Conclusions. The dead zone in our models possesses a mean magnetic field, significant Reynolds stresses and a steady local pressure 
maximum at the inner edge, where the outward migration of planetary embryos and the efficient trapping of solid material are possible. 

Key words. Accretion disks, magneto-hydrodynamics (MHD), instabilities, stars: planetary systems: formation, methods: numerical, 
turbulence 



1. Introduction 

Forming planets in a protoplanetary disk with a power-law sur- 
face density profile is difficult for several reasons. First, solid 
material on accumulating into meter-sized boulders quickly 
spirals to the s tar, transferring its o r bital angular momen- 
tum to the gas ( Weidenschillind 1977 : Nakagawa et alj T986t 



Takeuchi & Onl l2002t lYoudin & Chiand l2004t iBrauer et all 
20071) . Second, collisions between the constituents lead to dis- 



ruption rath er than growth when rather low speed thresholds 
are reached dBlum et alJ1998tlPoppe et al.ll999HBrum & Wurml 
120001) . Bodies in the meter size r ange are destroyed in impacts 
as slow as some cm/s dBenzl2000l) . Turbulence in t he gas readily 
yields collisions fast enough to terminate growth dBrauer et alj 
l2008ah . Third, Earth-mass protoplanetary cores are prone to ra- 
dial migration resulting from the tidal interaction with the gas, 
and in the classical type I migration picture quickly migrate 
all th e way to the host star dGoldreich & Tremaindl 1980c IWardl 
[T986h . 

All three of these problems could be solved by the pres- 
ence of local radial gas pressure maxima , which trap the drift- 
ing particles dHaghighipour & Bossll2003h . leading to locally en- 
hanced number densities and high ra t es of low-speed co llision 
dLvra et alJl2008L IBrauer et ai1l2008bt iKretke et al.ll2009l) . With 
sufficient local enhancement, one can envision the direct forma- 



tion of planetesimals via collapse under the self-gravity of the 
particle cloud, bypassing the size regime most susceptible to the 
radial drift. Furthermore, the radial migration of Earth-mass pro- 
toplanets can be slowed or stop ped by varying the s urface den- 
sity and temperature gradients dMasset et alj 120061) . Migration 
substantially slower than in the classical picture appears to be 
required to explain the observed exo planet population under t he 
sequential planet formation scenario dSchlaufman et alJl2009l) . 

The formation of local pressure maxima is governed by the 
radial transport of gas within the disk. The magne to-rotational 
instability or MRI fealbus & Hawlevlll99lL 1 19981) is currently 
the best studied candidate to drive such flows. Local shearing- 
box calculations show that the instability leads to long-lasting 
turbulence and to angular momentum transfer by magnetic 
force s, provided the magnetic fields are well-coupled to the 
gas dHawlev et alj 119951: iBrandenburg et all fl995l: ISano et al.l 
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Johansen et al I I 9h~ ideal MHD calculations 

have been perfor med in various astrophysica l contexts: pro- 
toplanetary disks dSteinacker & Hennindl200U lArlt & Riidigeil 
120011: lFromandl2005t iFromang & Nelson! I2006L l2009h . black 
hole accretion toruses (Hawlevl l2000l) . and galactic disks 
dDziourkevitch et al. 2004). All the global simulations included 
neither a physical magnetic diffusivity nor a physical viscos- 
ity. Meanwhile, from local shearing-box studies it is known 
that the strength of the saturated MR I turbulence depends crit- 
ically on the resistivity and viscosity dLesur & Longarettill2007l : 
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IFromang & Papaloizoul EzOO?: Fro mang et al.ll2007h . In particu- 
lar, whereas the molecular viscosity is small in protoste llar disks, 
the g a s is so weakly ion i zed in the cold disk interior dGammid 
119961: llgea & Glassgokl 119991: ISemenov et all 12004 that the 
Ohmic resistivity shuts down the linear MRI and prevents 
the dev e lopment of turbulence ( Sano et al. 1998t Sano & Stond 
l2002afllFle~ming & Stondl2003HTurner et al.ll2007h . 

In this paper we present the first global resistive MHD cal- 
culations to include the 'dead zone' where the rapid diffusion of 
the magnetic fields prevents magnetorotational turbulence. Local 
pressure maxima form in the calculations in two ways: at dead 
zone edges and in zonal flows. The dead zone edges yield long- 
lived rings of enhanced surface density near locations where a 
gradient in the ionization fractio n leads to a jump in the accre- 
tion stress riKretke&Linl 120071) . The zonal flows on the other 
hand result from local fluctuations in the Maxwell stress in the 
turbul ence, and lead to pres sure maxima with lifetimes of a few 
orbits dJohansen et al.ll2009l) . 

In the next section we describe our disk model and the choice 
of magnetic resistivity profiles. The third section presents the re- 
sults, starting from the global properties of the magnetic field, 
followed by the effects of the resistivity jump on the surface 
density. In the fourth section we discuss the interaction of the 
magnetic fields with the density rings and with radial minima 
appearing in the turbulent activity. Our main results are summa- 
rized in section 5. 




z ( c s /Q ) 

Fig. 1. Vertical profiles of magnetic diffusivity. Black lines 
show the profiles adopted for simulations of the dead zone, 
with ?7o = 2 • 10~ 6 (solid), r] = 7 • 10~ 4 (dot-dashed) and 
770 = 0.016 (dashed) (see Eq. [6j. Blue lines show the estima- 
tions of magnetic diffusivity made for the disk at 4.5 AU us- 
ing the simple gas-phase rea ction set (lOppenheimer & Dalgarnol 
ll974Hllgner & Nelsonl2006h together with dust grains. The solid 
line is for no grains, dotted for 0.1/j,m, dashed for 1/im and dot- 
dashed for 10/im dust grains. 



scribe the dead zone in the protoplanetary disk, we include the 
Ohmic dissipation in our models. The Ohmic dissipation is the 
largest non-ideal term in the in duction equation under typical 
dusty conditions dWardld 120071) . Ambipolar diffusion and Hall 
effect are not considered for the sake of simplicity. The vertical 
profile of magnetic diffusivity (see section 2.1) is adopted from 
separate chemistry calcullations and is fixed in space and time. 

We solve the set of MHD equations using 3D global simula- 
tions on a spherical grid (?•, 6, (f), 



dB 
~dt 



= V x [u x B - ?y(0)V x B], 



at p 



4irp 



[V x B] x B, 



dp 

at 



+ V • (pu) = 0. 



(1) 



(2) 



(3) 



The notation is the usual one. ^ is a point-mass gravitational 
potential and y/GMl is set to unity in the code units. We use 
the locally isothermal approach to describe vertically stratified 
disks, adopting P = c^(r)p(r,Q) with a sound speed c s = 
co/(r ■ sin(0)) and H/R = 0.07. The initial field setup is ex- 
actly the equilibrium solution, 



4(1 + a) 



1 



P = Po 



(r sin 0) a 



sin 



exp 



cos 6 
2c 2 sine 2 



(4) 
(5) 



where a — 3/2. 

We consider the inner part of the protoplanetary disk, which 
is heavier and warmer compared to a minimum solar-nebula. 
In order to mimic a 'dead zone' and the ionization thresholds, 
we adopt fixed magnetic diffusivity profiles 77(8). We use the 
estimations from dynamical di sk chemistry simul ations for the 
adopted disk parameters S, T dTurner et al.ll2007h . The surface 
density is S = (i?/3.75AU)~ 1/2 • 1700 g/cm 2 , and the temper- 
ature is T = (i?/3.75AU)" 1 • 280K, where R is cylinder radius. 
The units are [u$\ = 2ir AU/yr = 29.8 km/s for velocity and 
[77] = 2tt AU 2 /yr = 4.47 • 10 19 cm 2 /s for magnetic diffusivity. 
The models are listed in Table Q] 

Our models include the disk part from 2 to 10 AU. A purely 
azimuthal magnetic field is chosen as seed field for MRI tur- 
bulence, which is -P g as/-fmag = 25 everywhere in the disk. 
The Alfven limiter of CH m = 14co is applied. We use reflec- 
tive radial boundary condition, with buffer zones applied at radii 
2AU < r < 2.5AU and 9.5AU < r < 10AU. We apply the 
magnetic diffusivity within the buffer zones: ?7buffcr is 10 -5 at 
the radial boundary and decreases linearly towards the physical 
domain. The buffer zones have both damping of radial veloc- 
ity towards zero at the border, and diffusing away the magnetic 
eddies approaching the radial boundary. 

Periodic boundary conditions are applied both for azimuthal 
and vertical (i.e. 0) domain borders. 



2. Model Description 

The initial setup for the mo dels is very similar to those studied by 
IFromang & Nelsonl d2006l) for the ideal MHD case. In our mod- 
els, MRI-driven turbulence is operating in locally-isothermal 
disk, with a fixed spatial distribution for the temperature. To de- 



2. 1. Ionization Thresholds and Influence of Dust Grains 

An estimate of the magnetic diffusivity vs. height at 4.5 AU is 
shown in Fig. Q] The midplane diffusivity with dust grains ap- 
pears to be substantially higher then it is possible to include in 
the MHD simulations. The four blue curves from top to bot- 
tom are demonstrating the magnetic diffusivity in code units 
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Table 1. Model properties and midplane a-stresses inside (cka) and outside (cud) of the ionization threshold radius rth 
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Model 


Resolution 


rthlAU] 


VA 




no 






«D 


Time [years] 


Mm 


[256:128:64] 


4.5 


2.10' 


6 


7.10" 4 


1.9 ■ 10" 


3 


1.0 • 10~ 3 


880 


M m /2 


[128:64:32] 


4.5 


2.10" 


6 


7.10" 4 


1.3 ■ 10" 


-3 


2.1 • 10" 4 


960 


Mid 


[256:128:64] 


4.5 


2.10" 


6 


0.016 


1.6 ■ 10" 


-3 


2.9 • 10~ 4 


650 


Mi R6 


[256:128:64] 


6.5 


2.10" 


6 


7.10" 4 


2.9 ■ 10" 


-3 


1.6 • 10~ 3 


1100 


M idcal 


[256:128:64] 


none 


2.10" 


6 


2.10" 6 


5.4 ■ 10" 


-3 


5.4 • 10" 3 


682 * 



for the gas and dust grains of 0.1, 1 and 10 microns, and 
no grains. We have used the sim ple gas-phase reaction set of 
lOppenheimer & Dalgarnol ( 1 974]) tog ether with the grain surface 
chemistry of lllgner & Nelson ( 20061) for the classical dust to gas 
ratio. Ionization by stellar X-rays, cosmic rays and long-lived 
radionuclides is included. The penetration depths are assumed 
8g/cm 2 for the X-rays and 96 for the cosmic rays. 

The exact calculations of chemistry and dust behavior in the 
thermally evolving global disk is a hard task with many free pa- 
rameters. After planetesimals form, the dust mass fraction will 
be lower than the interstellar value. We shall b ear in mind that 
CRP stopping depth can be as low as 36g/cm 2 dGlassgold et al] 
l2009l 

Here we simplify the situation and adopt the following time- 
independent vertical profile of magnetic diffusivity, 



r\ = ?7o exp 



l6-1 



1.55 



(6) 



where 770 is the midplane value of magnetic diffusivity. In Fig.Q] 
black lines show the profiles adopted for our simulations, with 
??o = 2 • 10" 6 (solid), 770 = 7 • KT 4 (dot-dashed) a nd 770 = 0.016 
(dashe d ). To simulate the s ituation mentioned in iRretke & Linl 
J2007I) : iKretke et al] (120081) . we drop the radial dependence of 
j] and suggest that r/o makes the jump at the chosen threshold 
radius r t h- For example, in model Mir (section 2.2, Table 1), 
the Ohmic diffusion is following a black solid line in Fig. Q]for 
radii from 2 to 4.5 AU, and a dot-dashed line for radii from 4.5 to 
10 AU. The magnetic diffusivity in the disk without dust grains 
is low, so that we have MRI-active region from 2 to 4.5 AU. The 
dead zone begins where the dust grains are present. The exact 
location of the inner edge of the dead zone (i.e. r t h) depends on 
the properties of the star and the dust in the disk, and can vary 
from object to object. Our choice for locating the threshold r t h is 
quite arbitrary. We place r t h roughly in the middle of our inner 
global disk patch (TableQ}. 

2.2. Set of Simulations 

Models in Table Q] have an ionization threshold posed either at 
4.5AU or 6.5AU. Midplane values for magnetic diffusivity are 
noted in Table Q] as tja and ?/d (Active' and 'Dead') for gas 
states inside and outside the threshold radius. The time duration 
of each model is given in years, and the mark * is given when the 
steady-state has not been reached. Vertical profiles for magnetic 
diffusivity follow Eq.|6]and are demonstrated in Fig.Q]with black 
lines. Each model combines two diffusivity profiles, except the 
run Mideal- 

In Table 1, notations are T for quasi-ideal MHD state with 
77 = 2 • 1CT 6 (Fig. [T] solid black line), 'R' for the gas disk 
with lOyum-sized dust grains (770 = 7 • 10~ 4 , dot-dashed black 
line), 'D' for the case of l/urn grains (770 = 0.016, dashed black 
line). Our adopted magnetic diffusivity profiles for the disk with 
1/jm and 10/im-sized dust grains will allow the turbulent MRI 



layers beyond 3H and 2H, correspondingly. The peak values of 
blue curves in Fig.Q]are leading to unacceptable short time steps 
in resistive MHD simulations. We have observed that it is not 
convenient to compute the regions with 77 > 0.1AU 2 /yr with 
standard MHD codes, because of the dramatic shortening of the 
time step. For this numerical reason, we take the magnetic dif- 
fusivity slightly different as the chemistry models predict. Our 
adopted profiles of magnetic diffusivity (black curves, Fig.Q]) al- 
low to match the values of chemical models at 2 AU and remain 
above the numerical dissipation for region between 2H and 3H. 
Reducing of the magnetic diffusivity in the dead zone may in- 
fluence how fast the global magnetic fields are diffused into the 
dead zone, whereas the MRI modes are damped all the same. 

2.3. Calculation of Turbulent Stresses 

An important outcome of our simulations is the magnitude of the 
Reynolds and Maxwell st resses. To calculate t he lat ter, we use 
the approach described in iFromang & Nelsonl d2006l) for curvi- 
linear coordinates. The turbulent viscosity can be described as 
v = etc 2 / f2, where the main component of the a stress tensor is 



Tm + Tr 
IP) ' 



(7) 



or a = qr + ctM- The Reynolds and Maxwell stresses are cal- 
culated as 



T M = -((B'^B' r )/in), 

The mean pressure for azimuthal domain A</> is 
c s (r) 2 



(8) 



(9) 



(P(r))=c s (r) 2 £(r) 



3. Results 



A(j) 



pr sin (Q)ded(j). (10) 



Ae JAtj, 



In this section we describe our results and focus on two main is- 
sues. First, we study the time evolution and radial dependence of 
magnetic fields in our models. Secondly we study the formation 
of long-lived density rings which may or may not be able to trap 
solids in the disk and thus trigger the onset of planet formation. 
Table Q]represents the set of models. In 3.1 we discuss the issue 
of resolution. In 3.2 we describe the properties of global models 
with the emphasis on the evolution of the magnetic fields. In 3.3, 
the radial behavior of resulting Maxwell and Reynolds stresses 
is described. In 3.4 we explore the evolution of the pressure rings 
in time. In 3.5 we demonstrate the change of rotation and the tur- 
bulent properties of the gas in the rings and in the pressure bump 
at the inner edge of the dead zone. 
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Fig. 2. Inverse plasma (3. Top: model Mir at t — 300 years; Below: Mir/2 (halved space resolution) at the same time. Left: The 
change of inverse plasma /3 vertical distribution from initial (red) to the convex shape (black lines). The solid line stands for inverse 
plasma beta averaged within active zone (2.5 — > 4.5 AU). Dashed line is for the 1//3 averaged over the patch of the dead zone 
(5 — > 7AU). Vertical bars mark the disk height where the Elsasser number drops below unity (blue lines, solid for active and dashed 
for dead zones). Green lines border the midplane region with X c /A(j> < 8 (Eq.fTTI). Right: Radial dependence of vertically averaged 
1 //3 (solid line), 1//3 at the midplane (dotted line) and at z — 2H (dashed line). 



3. 1. Azimuthal MR I and the Issue of Resolution 

The MRI from a purely azimuthal magnetic field (AMRI) leads 
to non-axisymmetric perturbations. The radial displacements of 
the initial azimuthal field are enhanced due to the differential ro- 
tation. This leads to the appearance of field components B r and 
turbulent B$. The excess of magnetic pressure and the buoyancy 
lead to the generation of the vertical magnetic field component. 
The linear analysis has been done in Balbus & Hawley (1992). 
The critical wavelength for AMRI in units of the azimuthal grid 



X c /A(j) = 2tt 



16 2 

15/3 



c /A</>, 



(ID 



which follows from Eq. (15) in lHawlev et al.l d 1995b . When 
k ■ Va ~ and | k/k z | is in the right range, then the growth 
rate of the non-axisymmetric modes is greatest. In contrast to 
the magneto-rotational instability of vertical magnetic field, the 
vertical wavenumber by AMRI is not constant and increase with 



time during the linear stage of MRI. The largest total field am- 
plification is expected for the modes with largest possible \k z \. 
As a consequence, in numerical simulations of AMRI it may 
be impossible to resolve all growing wave -numbers and the to - 
tal amplification is limited by the gri d size (lHawlev et alj|1995l) . 
Nevertheless, the numerical study in lHawlev et al. ( 19951) shows 
that for effective resolution above 8 grids per critical wavelength 
in azimuthal direction the saturation of magnetic energy is only 
weakly affected by the resolution. Note, that the resolving of az- 
imuthal critical wavelength is important. In Fromang & Nelson 
(2006), the resolution of more then 5 grids per wavelength has 
been suggested as sufficient. 

All our runs are made with the initial uniform plasma beta 
of 25. Following Eq. [TTJ we have X c /A<p = 10.47 every- 
where in the MRI-active disk for the models with resolution 
of [256 : 128 : 64]. Model M m / 2 with halved resolution has 
A c / Acf> = 5.24. Note, that the Ohmic dissipation poses an addi- 
tional limitation for the excited MRI wavelength. 



Dzyurkevich et al.: Trapping Solids in 3-D Global MHD Simulations 



5 



Model M_IR, t = 300 years 



.2f 



7min(A) = 5.75lb-IC 
: min(A) = 0.0159 
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Model MJR6, t=300 years 




Radius TAUl 



Model M-TDEAL, t=300 years 




3 4 6 B 

Radius [All] 



Fig. 3. Elsasser number A = v\ z /rfl for models Mir (top), 
Mir6 (middle) and Mir (bottom). The Elsasser number is plot- 
ted in logarithmic color scale for edge-on view of the disk. 
Yellow-green color marks weak Bq magnetic fields which are 
stable to MRJ. Blue is marking the regions with the vertical mag- 
netic field sufficiently strong to launch the MRJ. The transition 
to stability at v\ z /rfl = 1 is indicated with a black line. Slices 
of turbulent magnetic fields are taken for t = 300 years (orbits 
at 1 AU) for each model. 



The effective resolution of 10.47 is holding during the linear 
stage of AMRI. The azimuthal magnetic field is breaking into 
filaments of opposite signs with Ae < A r < < A^ in MRI-active 
zones. The inverse plasma beta grows to 100 at the midplane 
and the decreases below unity in the upper disk layers. The ef- 
fect of expelling the magnetic field into the corona becomes vis- 
ible after 30 local orbits (section 3.2, Fig.[6]l. Note, that in low- 
resolution tests with A c /A0 = 2.62 (Mideal with resolution 
of r : : <j) = [64 : 16 : 8], excluded from Table 1), there is 



no MRI exited and the azimuthal magnetic field remains in its 
initial shape for few hundred years. 

Fig. [2] demonstrates how the inverse plasma beta, 1/(3 = 
Pmag/Pgaa, changes due to MRI from P mag /-Pg as = 1/25 (red 
dotted line) to the convex shape (black lines). The solid line 
stands for inverse plasma beta averaged within the active zone 
(2.5 — > 4.5 AU). The dashed line stands for the 1/(3 averaged 
over the patch of the dead zone (5 — > 7AU). In the midplane 
we find the minimum of magnetic pressure, with P mag /Pg as = 
0.01. The plasma beta is reaching 1 at 2.8 H both in model Mir 
and in the low-resolution run Mir/2 (Fig- 12 left). The resulting 
vertical profile of the magnetic pressure is very similar to those 
shown in Fromang & Nelson (2006). It is remarkable, that the 
dead zone builds up the same vertical distribution of magnetic 
pressure as the active zone, predominantly due to the smooth 
azimuthal magnetic field component. Radial dependence of the 
inverse plasma beta (Fig. [2] right) in the normal resolution run 
Mir shows that upper layers possess the constant P ma g/Pg a s, 
whereas the midplane layers, i.e. from midplane to 2H, are os- 
cillating and slightly decrease towards the inner radius within 
the active zone. The dotted line in Fig.|2](top right, model Mir) 
shows that the magnetic pressure falls to zero at r = 5.5AU at 
the midplane. The reason is a diffusion of the mean azimuthal 
magnetic field from the active zone into the dead zone. This dif- 
fused field has the opposite sign to the primordial field. Its time 
propagation into the dead zone is shown in Fig. |4] (section 3.2). 
The low-resolution model Mir/2 shows that the expelling of the 
azimuthal field into the corona is not reaching the same extend 
as in the normal resolution model for radii r > 6AU. 

This vertical re-distribution of the azimuthal magnetic field 
affects the effective resolution. In the left panels of Fig. [2] we 
adopt solid lines for active and dashed lines for dead zone val- 
ues. Green vertical bars in Fig.|2]mark the midplane region with 
A c /A0 < 8 (Eq. rTn >. The blue vertical bars show the disk 
height where Elsasser number A = v\ z /r)£l drops below unity. 
T he criterion for MRI i nstability A > 1 has been introduced 
in ISano & Stonel d2001l) for the case of non-ideal MHD with 
Ohmic dissipation. After 300 years, the B$ vertical profile is 
changed so much that the midplane layers are resolved only with 
A c / > 5, for example in the active zone (2.5 AU to 4.5 AU) in 
normal resolution model (green bars, Fig. |2j. On the other side, 
model Mir/2 becomes well-resolved in the layers \z/H\ > 2. 
Effective resolution of X c /A<f> > 16 is reached in the active lay- 
ers above the dead zone, where vertical MRI is launched outside 
of the A = 1 line. Interesting to note, that the Elsasser num- 
ber A = v\ z /r)Q. drops below unity roughly at the same height, 
when we compare normal and low resolution models (Fig. [2] 
top left and bottom left). When looking for numerical values of 
A c /A0 get at the location of blue bars in the active zone, we 
find A c /Ac/> = 6 for Mir and 4 for Mi R / 2 . These numbers are 
only approximate values, because it is difficult to calculate them 
accurately at the height at which A = 1. All in all, there are 
surprisingly small differences between Mir and Mjr/2 models. 
The instability occurs in both cases and leads to similar physics, 
though the speed of the total field amplification is slower for 
Mjr/2 (Table 1). 

Fig.[3]demonstrates the Elsasser number and A = 1 criterion 
in the (r, 0) plane of azimuthally averaged disk. The solid black 
line of A = 1 marks the locations where the regeneration of 
vertical magnetic field cannot be efficient enough. Models Mir 
and Mir6 have their most active MRI layers at 2.5H above the 
'dead' midplane. Comparison of the contour plots of the Elsasser 
number for Mir, Mirb and Mideal shows, that MRI-active 
zones have very similar appearance. At the midplane of the ac- 
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Fig. 4. Temporal evolution of azimuthal magnetic field as a 
function of radius for Mir model. Azimuthally averaged 
field is demonstrated for horizontal slices atO = tt/2 + n ■ cq 
or at z = n ■ H, where n = 0, 2. In model Mir, the MRl-active 
zone reaches from 2.5 to 4.5 AU and layer z = is MRI-'dead' 
between 4.5AU and 10AU. The period of the B<f, sign reversals 
is about 150 years (30 local orbits). The black dashed line repre- 
sents the time of 10 local orbits for each radius. 
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Fig. 5. Temporal evolution of azimuthal magnetic field as a 
function of radius for Mir6 model. Azimuthally averaged 
field is demonstrated for horizontal slices at Q = tt/2 + n ■ cq 
or at z = n ■ H, where n = 0, 2. In model Mir6, the MRI- 
active zone reaches from 2.5 to 6.5 AU. The period of the B^ 
sign reversals is about 150 years (30 local orbits), same as in 
Fig. [4] Model Mir6 demonstrates that sign reversal happens not 
only in time, but along the radius as well. The black dashed line 
represents the time of 10 local orbits for each radius. 



tive zone, there are yellow-orange areas of low Elsasser number 
at r - 4AU (M m , M m6 ), r ~ 4AU (Mideal), and r ~ 5.3AU 
(MiR6), which are corresponding to the enhanced gas pressure. 
In section 3.4 we describe the formation of the density rings at 
these locations in more details. Conditio n V^/rjfl > 10 is ful- 
filled in the whole disk in every model (iTurner & Sand 12008b 
and the magnetic fields may be pumped into the dead zone from 
the active layers. 

3.2. Oscillating Magnetic Fields and Saturation 

After most of the initial azimuthal magnetic flux has been shifted 
to the upper disk layers, the two scale heights to the adjacent 
midplane develop the oscillating axisymmetric magnetic field. 
For model Mir, the oscillations last for over 400 years and then 
decay. In the MRI-active zone, i.e. between 2.5 and 4.5 AU, the 



Btf, sign-switching occurs within about every 120 to 150 years. 
In Fig. |4] we demonstrate the time evolution of the azimuthally 
averaged as a function of the radial distance for the Mir 
model. The fluctuating part of the azimuthal magnetic field is 
roughly ten times stronger than the mean part, B^ ~ 10(B^). 

This property appears to be typical for MRI, as the ob serva- 
tions o f the magnetic fields in the galactic disks show (see lBeckl 
(l2000t) for review). When the active zone is stretched up to 6.5 
AU, it is not affecting the time period of the B$ sign reversals. 
Model Mir6 demonstrates that sign reversals occur not only in 
time, but also along the radius, as shown in Fig. [5] The first pos- 
itive B^-stripe is starting at 60 years at 3 AU and progresses to 
6 AU in 160 years (Fig. [5). The stripes of B^ in (r, t)-plane are 
stretched along the line of local orbit, which is plotted with black 
dashed line for £i OC al = 10 in Fig. [5] At later times, the mixing 
and interaction of the waves is breaking radial B<j, reversals into 
less regular oscillations both in Fig. [4] and in Fig. [5] The field 
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Fig. 6. Horizontally averaged magnetic fields components (B r ), (Bq) and (B^) as function of disk height and time. Blue and red 
colors are tracing the positive and negative fields, green is always zero. (B^) sign-reversals 'swim' from midplane to the disk corona 
('butterfly diagram'). Above: model Mir, averaging of the magnetic fields has been done within the MRI-active zone (between 2AU 
and 4AU). Below: model Mir, averaging of the magnetic fields has been done within the dead zone (between 6AU and 8AU). 



diffusion in the dead zone outside of rth = 4.5 AU (model Mir) 
and outside of r t h = 6.5 AU (model Mir6) does not follow the 
line of local orbit, oc r 3 / 2 , but propagates with diffusion time 
oc r 2 . For the Mir run, the region with r > 4.5 AU develops the 
negative azimuthal magnetic fields due to diffusion of the mag- 
netic field. After 600 years there is a weaker wave of positive 
B^, which has a shorter time period. 

Oscillations in the sign of B^ at the z > 2H appear less 
clearly, compared to the azimuthal magnetic field at the mid- 
plane (Fig. 3] Fig.[5]l. A comparison with Mid shows why. The 
reason is the interaction of MRI waves in the upper layers be- 
tween active and dead zones. For the case of a very thick dead 
zone (Mid), the sign reversals in B^ at z = 2H show most pro- 
nounced intervals. This is due to the fact, that MRI can be best 
exited between 1H and 2H layers of the active zone. In the model 
Mir, the layers above and below the dead zone are turbulent and 
interacting with MRI modes at same height in the active zone, 
what leads to a more irregular picture in oscillations of B ( p(r 1 1) 
at z > 2H. The color-coded presentation of B^ as a function of 
(z, t) reveals a butterfly diagram, if the azimuthal magnetic field 
is averaged within the disk region between 2 and 4 AU (Fig.|6l 
top). The low plasma (3 does not prevent the upper disk layers 
from being very turbulent, as one can see from left and middle 
panels for B r and Bq. Comparing the contours of dominating 
B,p(z, t) and other two turbulent field components B r (z, t) and 
Bq(z, t) demonstrates again that the MRI turbulence and verti- 
cal redistribution of azimuthal magnetic field are connected. The 



period of B^ oscillations is about 150 years, corresponding well 
to the radial changes of B^ sign shown in Fig. [4] and Fig. [5] 
When averaging over the whole radial extent of the disk, or at 
least within the dead zone, then the butterfly picture disappears 
(Fig. [6] bottom). 

Volume-averaging of the magnetic energy shows, that its to- 
tal value is oscillating in time, with a period correlated to the 
sign-switch in azimuthal magnetic field within ±2H relative to 
the midplane. Fig. Q shows the magnetic energy for each model 
in Table 1 and the corresponding total alpha stresses. The os- 
cillations of energies in time can be clearly correlated with the 
butterfly diagram. Oscillations are weakly visible in the total al- 
pha stress (Fig. |7J right). The magnetic energy curves reach a 
constant value in model Mir , which has the smallest active zone 
(from 2.5 to 4.5 AU). Model Mid looses the total magnetic en- 
ergy continously during 400 years due to higher magnetic dis- 
sipation. Model Mid reaches a steady-state when stresses and 
magnetic energy remain unchanged from 400 to over 600 years. 
The longer is the MRI-active domain, the longer it takes for the 
simulation to reach the steady-state. Models Mir6 and Mideal 
hold oscillatory (non-stationary) magnetic fields for 900 years. 

The closed boundary conditions enforce the conservation of 
the total flux in the domain. Fluxes of vertical magnetic field re- 
main zero through the whole simulation. The effect of the bound- 
ary choice on the butterfly diagram remains to be investigated in 
future work. The local box simulations, made for open vertical 
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Table 2. a-stresses inside (A) and outside (D) of the ionization threshold radius for four layers above the midplane 
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Fig. 7. Top: Total magnetic energy evolution; Bottom: Total al- 
pha stress evolution for models Mir, Mid, Mir6 and Mideal- 
The oscillations are correlated with the 'butterfly diagram' . 



bound aries, show the butterfly picture as well dTurner & Sand 
l2008h . 



3.3. Maxwell and Reynolds Stresses 

Radial inhomogeneity in turbulent viscosity has been suggested 
as the mechanism to produce the pressure maximum, which is 
efficient in dust trapping, and therefor important in the planet 
formation theory. In our models, the turbulent viscosity is driven 
by MRI and the inhomogeneity in turbulent stresses appears nat- 
urally as the result of simulations, when we include the sharp 
gas ionization threshold. Indeed, we find a density bump forming 
behind the ionization threshold in our simulations (section 3.5), 
and a corresponding jump in turbulent a stress. 



Fig. 8. Snap-shots of turbulent a stress in (r, 0) (top) and (r, (f) 
(bottom) slices for models Mir for data output at t = 610 years. 
Black line indicates the ionization threshold. 



The time evolution of the Maxwell stress Tm (r, t) (Eq. |9]l 
is demonstrated in Fig. [9] There is a weak Maxwell stress of 
about 10~ 5 in the dead zone, which can periodically become 
negative. One can see the sharp border in Tm(t, t) between the 
active and the dead zones in slices for z = 0, z = 1H and z = 
2H. The traces of sign reversals in the azimuthal magnetic field 
are also visible in the Maxwell stress. Fig.|9]shows that exactly 
at 150 years the Maxwell stress reaches its maximum of 10 _1 , 
when calculated in units of initial pressure (P(r)) (Eq. |9j. Later 
on, the saturation of MRI sets in and the total stress is between 
10 -3 and 10 -2 (see Table 1). The dark-orange filaments of very 
weak Maxwell stress in the active zone, ±10 -7 , correspond to 
the location where the reversals of axisymmetric azimuthal field 
happen. The weak Maxwell stress is located at r — 3.5 AU for 
many years, what appears as a systematic stripe when looking at 
z = and z = 1H horizontal slices of Tu(r, t) in Fig. [9] This is 
a location where the density ring is created (more in section 3.4) 
and also most of B$ reversals take place during the time period 
from 200 to 700 years. Negative values of Maxwell stress appear 
at 3H above the midplane. This is the region of low plasma beta, 
and the turbulence at this height is no more MRI-driven. 

In the dead zone, the value atotal ~ 10~ 3 is due to Reynolds 
stress contribution. In order to provide the understanding of how 
the turbulence there looks like, we present a snap-shot of turbu- 
lent a stress for the Mir model (Fig. [8j. The dead zone is filled 
with vertical pillars of the a stress of opposite sign. In the (r, 4>) 
plane, those pillars look like tightly-wrapped spirals. The spiral 
waves are launched from the dead-zone edge (r t h = 4.5 AU), 
where the non-axisymmetric fluctuations in all velocity com- 
ponents are MRI-generated. The weak spiral structures can be 
found in the gas density as well. 
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Fig. 9. Model Mir: temporal and radial evolution of Maxwell stress. Horizontal slices of Tm(t, t) are taken at0 = 7r/2 + n ■ Co 
or z = n ■ H, where n = 0, 1, 2, 3. The color scale shows the azimuthally averaged Maxwell stress in local gas pressure units. 
Dark-orange filaments (±10 -7 ) in the active zone correspond to the location where the reversals of axisymmetric azimuthal field 
happen. Weak Maxwell stress of ±10~ 5 also occurs in the dead zone. 



We have calculated the turbulent stresses as the function of 
radius. The vertical averaging has been done separately for four 
midplane-symmetric layers: |0 — 7r/2| < cq, cq < |0 — tt/2\ < 
2c , 2c < |0 - ir/2\ < 3c and 3c < |6 - ir/2\ < 4c , 
where cq = Hj R = 0.07. The resulting a stresses for these lay- 
ers are given in Table [2] as ±1H, 2H, 3H, AH correspondingly 
and plotted in Fig. [10] and Fig. [TT] with solid, dotted, dashed and 
dot-dashed black lines. In Tabl^2] the mark * indicates that the 
steady-state has not yet been reached. Solid blue lines in Fig.fTOl 
and Fig. QT| represent the a(r) integrated over the whole disk 
thickness. The strongest total stress alpha of all models is ob- 
tained in model Mideal, which remains of the same order of 
magnitude with radius. 



The turbulent stresses have oc r~ 2 slope during the 'but- 
terfly' evolution stage in models MiR6 an d Mir. For the Mir 
model, the time averaging of the stresses between 250 and 600 
years results in a a oc r -2 . Afterwards, the turbulence in Mir 
reaches a 'butterfly' -free state (t > 600 years). When the tem- 
poral averaging is made between 400 and 800 years, instead of 
from 250 to 600 years, we obtain a\ which is constant with 
radius (Fig. ITTb . In the Mir model, the Maxwell stress is the 
main contribution to atotai only for layers 3H, AH outside the 
threshold r th = 4.5AU (Fig.HD. 

The gas density perturbations look like spiral waves, and 
the Maxwell stress is strongly reduced outside of the ionization 
threshold. Model Mid shows a steeper fall in Reynolds stress 
and in total a in the dead zone, compared to the Mir model. 
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Fig. 10. Reynolds, Maxwell and total a stresses for models Mideal (left) and Mir6 (right). Vertical averaging of stresses is made 
separately for four midplane-symmetric layers, ± lif (solid line), 2H (dotted line), 3iJ(dashed line), AH (dot — dashed line). 
Solid blue line shows the classical total stress. All stresses are normalized to the pressure at t = 0. Time average is made between 
850 and 1050 years for Mtr,6, and between 610 and 681 years for Mideal- 



The active layers above and below the dead zone are very thin 
and only marginally unstable to MRI. Model M ro has 1H, 2H 
and partly 3H layers which are deactivated in the dead zone. 
The pumping of the waves happens mostly from the MRl-active 
zone, and not from both active zone and adjacent layers as in 
Mir. Comparing models Mir and Mid, we conclude that the 
pumping of the hydro-dynamical waves into the dead zone is 
coming both from 'sandwich' active layers and from the active 
zone within r t h- 

The summary of turbulent a stress values for each layer 
above the midplane is presented in Table [2] The radial aver- 
aging has been done separately for active (A) and dead (D) 
zones. The turbulent a(A) stress is increasing from midplane to 
2co < 1 6 — 7r/2| < 3co, and drops in the fourth layer. The most 
prominent decrease in a stress from active to dead zone hap- 
pens within adjacent to the midplane |0 — ir/2\ < 2cq layers 
for models Mir and Mirq. In the case of very thick dead zone 



(Mid), the decrease in turbulent stress at the dead zone edge is 
significant in all three layers, |0 — tt/2\ < 3cq. 

3.4. Development of Pressure Maxima and Trapping of 
Solids 

Due to active secretion through MRI-active zone, the pressure 
bump at the inner edge of the dead zone is formed within a hun- 
dred of inner orbits. 

The change in the radial density gradient appears in the disk 
layers starting from the midplane and up to 3H. The piling up of 
the density behind the r t h is accompanied by a broad gap before 
the threshold. In addition, we find the rings of enhanced density 
within the MRI-active zone. The pressure bump at the inner edge 
of the dead zone and the rings of enhanced density within the 
MRI-active zone are formed due to different processes, which 
we discuss in section 4. Fig.[T2l(top) shows the changes of sur- 
face density for models Mir and MiRg. In our simulations, the 




Fig. 11. Reynolds, Maxwell and total a stresses for models Mir (left) and Mid (right). Vertical averaging of stresses is made 
separately for four midplane-symmetric layers, ± lif (solid line), 2H (dotted line), 3iJ(dashed line), AH (dot — dashed line). 
Solid blue line shows the classical total stress. All stresses are normalized to the pressure at t = 0. Time averages has been done 
starting from time of 400 years until the last data output. 



number of density rings depends on the extension of the active 
zone. In the model without a dead zone, there are three rings of 
enhanced density appearing within the domain before the quasi 
steady-state is reached (Fig. [12] bottom right). The thicker dead 
zone (model A/id) seems to be more efficient in piling up the 
higher bump: the maximum of the surface density peak is larger, 
when compared to the snap-shots of surface density in Mir at 
the same time. 

First, we consider the accumulation of the density at the in- 
ner edge of the dead zone. The pressure bump is fixed in time 
at the location behind the ionization threshold rth> as in Fig. IT3l 
(left). There are three stages in the evolution of the pressure max- 
imum at the inner edge of the dead zone (Fig. Q~3] middle), for 
example when considering model Mir. One can recognize the 
period of very fast mass excavation for the time from t = 
to 150 years. During t = 150 — > 600 years the peak of sur- 
face density is still growing at a roughly ten times slower speed. 
After t = 600 years there is no further increase of surface den- 



sity, what corresponds to the steady-state in the presence of sat- 
urated MRI-turbulence. The most straight-forward explanation 
for three stages in density excavation gives a time-dependent 
evolution of the Maxwell stress, since it governs the accretion in 
the MRl-active zone and in the active layers. The quasi-steady 
state is reached for models Mir and Mid after 600 years, and 
the maximum of the pressure bump at the ionization thresh- 
old and the minimum of the surface density in the gap remain 
unchanged. The density rings in the active zone appear to be 
long-living but less stable features. We observe the merging of 
two rings of enhanced density after 640 years in model Mir6 
(Fig. [T3l bottom). 

The radial pressure gradient is negative in the smooth un- 
perturbed disk, what leads to sub -Keplerian gas rota tion. Dust 
grains undergo an orbital decay (lAdachi et al.l 1 1976h . because 
they experience the 'head' wind. For example, the meter-size 
particle will migrate from 1 AU into the Sun within few hundred 
years (MMSN model). When local positive exponent of the disk 
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Fig. 12. Surface density after 600 years of evolutions. Red solid line stands for the initial surface density profile, black solid line 
represents the final surface density. Red dotted line indicates the cosmic ray adsorbtion depth 100g/cm 2 . Blue lines show the 
contributions of every scale height in the surface density. Vertical dashed line shows the inner edge of the dead zone. 



midplane pressure appears, the dust grains may experience 'tail' 
wind and the hydrodynamical d rag will lead to their outward 
migration (N akagawa et al.|[l986l) . The criterion for outward mi- 
gration of the dust grain due to the gas drag is 

p>-q/2 + 3/2, (12) 

where p = d\og £ gas /<ilogr and q = 2d log c s /d log r. Note 
that the factor 3/2 is the normalized shear. The p lanetary em- 
bryos can be stopped in such pressure traps as well dZhang et al.l 
120081) . Criterion for ou tward migration o f the protoplanetary 
cores has been given in llda & Lid (l2008h and is based on the 
migration rate of planets dTanaka et al H2002I) . When curvature 
effects on the Lindblad resonances have been included, the con- 
dition for outward migration for embryos is 

p > -0.80(7 + 2.52. (13) 

The right panels in Fig. [l3]show where criteria for outward mi- 
gration are fulfilled in our models. Black lines stand for criterion 
given in Eq. Q~2] and the red dashed line is for Eq. [T3] The out- 
ward migration of planetesimals is possible within the pressure 
bump at the inner edge of the dead zone (Fig.Qj). Conditions for 
embryos are more difficult to satisfy. Not all density rings within 
the active zone provide sufficiently strong positive pressure gra- 
dients, so that planetary embryos cannot be stopped there. The 
exception is model Mir . There we obtain two pressure traps at 
r = 3 AU and at r = 4.5AU, where large bodies can be retained. 



3.5. Super-Keplerian Rotation and Similarity to Zonal 
Flows 

When a quasi-steady state is reached, all time derivatives can 
be neglected and the Navier-Stokes equation for radial velocity 
gives: 



dr 



1 dc 2 s p 



P 
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lorcntz ■ 



(14) 



In our locally-isothermal simulations, the temperature is con- 
stant on cylinders. The steady-state solution is then given as 
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It follows from the last equations, the gas may reach purely 
Keplerian rotation when the density profile is locally changed 
to p oc r and the Lorentz forces are negligibly weak. 

The vertical profiles of turbulent velocity dispersion and the 
mean radial drift velocity AU = u$ — Ukcpi are plotted in 
Fig. [14] (Mir) and Fig. Q3] (Mire). The root-mean-square tur- 
bulent velocities can be described with the same shape of the 
vertical profile both in the active and in the dead zone. The tur- 
bulent velocity dispersion reaches 0.4 in the corona. This agrees 
with the result of Fromang & Nelson (2006). The top left pan- 
els in Figs. [14] and Fig. [15] show the root-mean-square turbulent 
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Fig. 13. Relative surface density S/£n, for models Mir (above), Mid (middle) and Mirq (below). Left: Colors show the evolu- 
tion of density bumps (yellow) and gaps (green). Middle: time development of max(E/I]o) of the bump at r = [3 : 6]AU and 
min(S/So) of the gap in r = [2.5 : 5]AU. Right: Criteria for outward migration (Eqs. fl2l [T3l l. Rings of enhanced density are 
formed in the active zone and the outmost density ring is located at the ionization threshold, marked as a vertical line in both left 
and right panels. 



velocities averaged separately for active (A) and for dead (D) 
zones. It is interesting to note, that the levels of yj < u^/c 2 >, 
^/< u 2 /c 2 > at the midplane are not very different for the two 
zones. Midplane values vary from 0.03 for \/ < Uq/c 2 s > to 0. 16 
for < u^/c 2 >. As expected, the perturbations of rotational 
profile are higher in the active zone. The radial dependence of 
\J < u@/c 2 >, y/ < u 2 . /c 2 > shows that the turbulent dispersion 
is significantly reduced within the pressure bump at the inner 
edge of the dead zone (Fig. [14] and Fig. [15] bottom right). 



The time-average of AU shows that the super-Keplerian ro- 
tation at the density rings is a long-lasting effect. Relatively large 
dust particles (i.e. for Stokes num ber St = 1) will m igrate out- 
wards with the velocity v r = AU (Kla hr & Linl200lT) . The areas 
of outward migration are up to 0.5 AU broad and outward veloc- 
ities reach 0.15 km/s (Fig.[T4"land Fig. [15] top right). 

There are certain similarities between the de nsity bump and 
rings w e found and the zonal flows described in ljohansen et al.l 
(120091) . The enhancement in density has been reaching 10 pres- 
sure scale heigh ts both in large loc al-box simulations and in the 
global model of iLvra et al.1 d2008l) . We estimate the rings to be 
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Fig. 14. Radial and vertical turbulent properties in model Mir. Top left: root-mean-square turbulent velocities are averaged sepa- 
rately for active (A) and for dead (D) zones. Top right: Outward migration velocity. Panels on bottom right and bottom left show 
time-averaged density, Maxwell stress, magnetic energy and r.m.s. velocities at the midplane. Time average is from 600 to 800 years 
(steady-state stage in model Mir). 



roughly 1.5 AU broad, what gives maximum 6 pressure scale 
heights. We show the radial correlations of the following pa- 
rameters: AU = it0 — (7Kcpi for various heights in the disk, 
Sp/po = (p — po)/po, magnetic pressure B 2 /B 2 , and Maxwell 
stress BrB^/Pinit (Fig. [14] Fig.[l5j bottom left). The gas around 
the rings of enhanced density within the active zone in MiRg 
and Mir models shows similar turbulent propert ies to the zonal 
flows and corresponding density maximum as in lJohansen et alJ 
(l2009h : The magnetic pressure and Maxwell stress are strongest 
where the density minima are excavated, the maximal AU is 
half -phase shifted. The maximal deviatio ns from Kepler velocit y 
have same amplitude as in zonal flows in lJohansen et al.l d2009t) . 
but the relative amplitude of density Sp/ po = (p — Po)/ Po is ten 
times larger in our models. The maxima of magnetic pressure 
B 2 /Bq and Maxwell stress B r B^ / Pj n i t are located between the 
rings and are about te n times stronger compared to the values in 
iJohansen etaD (l2009h . 



In addition, the pressure bump at the ionization thresholds 
has significantly less turbulence; the Uq/c 2 has a value only 
about 0.05 and not 0.5 as outside of the density bump. In the 
case of Mir, the density ring in the active zone is a 'calm place' 
as well. Model Mtr6 keeps the butterfly pattern of the azimuthal 
fields much longer. This is the reason why the turbulent veloc- 
ities decrease in the density rings weaker in the case of MiR6, 
compared to Mir. 

4. Synthesis: Connection Between Pressure 
Maxima and 'Butterfly' Structures 

4. 1. Density Rings in the MR I- Active Region 

In orde r to provide the comparis on with previous global simu- 
lations (iFromang & Nelso"nll2006l) . we have done the fully MRI- 
active disk model. The turbulent and magnetic properties of the 
fully MRI-a ctive disk model Mideal are very close to those 
presented in lFromang & Nelson! d2006l) for run S4. The toroidal 
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Fig. 15. Radial and vertical turbulent properties in model Mir,6. Top left: root-mean-square turbulent velocities are averaged sepa- 
rately for active (A) and for dead (D) zones. Top right: Outward migration velocity. Panels on bottom right and bottom left show 
time-averaged density, Maxwell stress, magnetic energy and r.m.s. velocities at the midplane. Time average is from 900 to 1100 
years. 



magnetic field is expelled from the midplane into the upper 
disk layers within the linear stage of MRI turbulence. The peak 
value of the volume-integrated turbulent a stress is reaching 
0.019 in our model Mi DE a l (Fig. 01 and 0.013 in the S4 model 
(iFromang & Nelso"nll2006l) . At the end of the simulation, the tur- 
bulent a stress is 0.005 in gas pressure units. 

In our simulations, the MRI turbulence evolves through three 
stages: (a) linear growth, (b) oscillatory saturation regime and 
(c) non-oscillating steady state. The stage (b) is best to observe 
if the dead zone is included. The oscillations are regular and can 
be registered in total magnetic energy and in turbulent a stress 
time evolution. The dominating azimuthal magnetic field com- 
ponent in the MRI-active zone switches its sign with a period 
of about 150 years or 30 local orbits, and within the radial ex- 
tent of 1.5 AU. The sign reversals of the azimuthal magnetic 
field with respect to the midplane, known as butterfly diagram, 
have also been observed in local shearing-box simulatio ns of the 
stratified disk dTurne r et al. 2 007t Hohansen et ai]l2009h . In case 



of the global simulation (model Mideal), averaging over the 
whole eight AU of the MRI-active domain brings a very irregular 
butterfly picture. Numerous reversals of the azimuthal magnetic 
field along the radial extent lead to seemingly irregular peaks in 
a(r) and in the magnetic energy. 

We observe in our models, that the reversals in B$ come 
along with the density rings. Thus, it is important to understand 
what causes radial and temporal reversals in the magnetic field, 
and what determines a period. We observe that the amplitude of 
the oscillations in magnetic energy is higher, if we increase the 
radial extent of the active zone (models Mir and MiRg). In ad- 
dition, it is important to know how long the 'butterfly' structure 
can survive. The magnetic field is oscillating over the whole du- 
ration of the local-box simulations. In our global runs, the life- 
time of the oscillating stage (b) seems to depend on the extent 
of the MRI-active zone (Fig. 0. The thickness of the dead zone 
also influences the life-time of the oscillatory regime (b), as we 
found from a comparison of magnetic energy curves for models 
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Mir and Mid. The MRI waves in active layers above and below 
the dead zone are interacting with the turbulent magnetic field 
inside the threshold radius. If the layers above and below the 
dead zone are as thin as in our model Mid, the turbulent a stress 
in the active zone appears to be lower, the same applies to the to- 
tal magnetic energy. The oscillations of the azimuthal magnetic 
field have ceased after 400 years, and after 600 years in the case 
of thinner dead zone (model Mir, Pig. EJ. The turbulent a stress 
is roughly constant with radius during the non-oscillating steady 
stage (c). 

The inner radial boundary with its resistive buffer (2AU < 
r < 2.5AU) and the inner edge of the dead zone r t h enclose the 
MRl-active zone in our models Mir, Mid and Mnte. We ob- 
serve a butterfly diagram in the MRI-active zone of our models, 
when plotting (B^(z)) oc t. There are no oscillations of (B,p(z)) 
in or around of the dead zone. We obtain the clearest oscillations 
of at the midplane of active zone, |9 — 7r/2| < cq. The B^- 
sign reversal happens every 150 years and within 2.5AU < r < 
4AU (models M m , M m6 ) and 4AU < r < 5.5AU (M m6 ), i.e. 
within every 1.5 AU. Those reversals are stretched along the line 
of the local orbit in (r, t) space. 

The rings of enhanced density appear already at t = 150 
years, roughly at the time when the linear AMRI breaks into a 
nonlinear regime and the 'butterfly' is initiated. During the os- 
cillatory stage (b) of the MRI evolution, there is a radial depen- 
dence in turbulent a stresses according to a(r) oc r~ 2 . At the 
inner radii, the a stress is high and it leads to the effect of fast 
local excavation of the density and accumulation of it at some 
outer radius. The radial reversals of the azimuthal magnetic field 
are aligned with the rings of enhanced density. The radial loca- 
tion of magnetic field reversals remains constant over hundreds 
of years. At the same location, we find the stripes of weakest 
Maxwell stress ~ 10~ 7 (Fig. 1 1). As soon as the ring of density 
is created at a certain location, there is a corresponding change in 
the rotation and in the shear. On the one hand, over-density leads 
to the local P mag /P gas relation which can be too low to excite 
AMRI. On the other hand, the change in the rotation reduces the 
shear and it leads to local stabilization of MRI within the den- 
sity ring. The consequence is that the rings of enhanced density 
are less turbulent compared to the density minima between the 
rings. In contrast to the local-box studies dJohansen et al.ll2009t 
iGuan et aLll2009t) . there is more then one density ring forming 
in the case of a longish MRI-active zone. The merging between 
rings is possible. In the Mideal model, we observe the appear- 
ance of three weak density rings. 

4.2. Pressure Maximum at the Inner Edge of the Dead 
Zone 

The formation of the pressure bump at the dead zone edge is 
not directly correlated with the 'butterfly' magnetic structures 
within the active zone. The strength of the turbulent viscosity in 
the MRI-active zone determines the time scale to form such a 
pressure bump. The extent of the MRI-active zone has an effect 
on the value of the total turbulent a stress in the active zone. 
As it follows from Table 1, the turbulence in models Mideal 
and MiR6 provide the largest stresses of 5.3 ■ 10 -3 and 2.9 • 
10" 3 , followed by M m with 1.9- 10~ 3 andM ro with 1.6- 10~ 3 . 
The radial extent of the active zone was reduced from 8AU to 4 
AU and to 2AU. This sequenc e is co nsistent with the results of 
IGuan et al ] (120091) : lBodo etal.1 (120081) . where a similar decrease 
of total alpha stress is demonstrated when the size of the local 
box is reduced. 



The jump in the turbulent viscosity is responsible for the for- 
mation of a pressure trap at the dead zone edge. This jump is 
strongest at the midplane. The intensity of MRI-driven turbu- 
lence grows strongly with the disk height. The locally calculated 
stresses are different for each pressure scale height layer: there 
is up to 1 order of magnitude difference when comparing the 
values for midplane and top (|0 — 7r/2| > 3co) layers of the 
active zone. As expected, the midplane value of Maxwell stress 
decreases drastically in the dead zone, but the Reynolds stress 
does not 'feel' the presence of the dead zone. We observe only a 
slight bump in the Reynolds stress at the location of the density 
maximum. Within the dead zone, there are spiral density waves 
propagating from the inner edge outwards. The vertical velocity 
dispersion is non-zero as well. We conclude that the resulting 
a of about 10 -3 is due to the waves pumped vertically from 
the active layers and radially from the active zone through the 
threshold radius. For the disk evolution, it is important to have 
a significant residual a in the dead zone in order to reach the 
quasi-steady state without getting unstable (i.e., a gravitational 
instability of the density ring). 



5. Conclusions 

We have presented the results of the first global non-ideal 3D 
MHD simulations of stratified protoplanetary disks. The domain 
spans the transition from the MRI-active region near the star to 
the dead zone at greater distances. The main results are as fol- 
lows. 

- As suggested in lKretke et al.l d2008l) . the height-averaged ac- 
cretion stress shows a smooth radial transition across the 
dead zone edge. The stress peaks well off the midplane at 
2co < 1 6 — 7r/2| < 3co- Consequently, averaging over the 
full disk thickness yields only a mild jump, despite the steep- 
ness of the midplane radial stress gradient. 

- Weak accretion flows within the dead zone are driven mainly 
by Reynolds stresses. Spiral density waves propagating hor- 
izontally produce non-zero velocities and angular momen- 
tum transport near the midplane, apparently with little asso- 
ciated mixing. The dead zone midplane vertical velocity dis- 
persion |w e | is 0.03 times the sound speed, two to three times 
less than the radial and azimuthal components (Fig. IT5Tl. The 
waves are pumped both by the active region near the star and 
by the active layers above and below the dead zone. In model 
Mid where the dead zone is very thick, the pumping from the 
surface layers is weak and the stress falls steeply with radius 
within the dead zone (Fig. [TTJi. 

- The excavation of gas from the active region during the linear 
growth and after the saturation of the MRI leads to the cre- 
ation of a steady local radial gas pressure maximum near the 
dead zone edge, and to the formation of dense rings within 
the active regi o n, res embling the zonal flows described in 
Uohansen et all (120091) (Fig. [12] Fig. [13). Super-Keplerian 
rotation is observed where the radial pressure gradient is 
positive. The corresponding outward radial drift speeds for 
bodies of unit Stokes number can exceed 10% of the sound 
speed. The pressure maxima are thus likely locations for the 
concentration of solid particles. 

- The turbulent velocity dispersion, magnetic pressure and 
Maxwell stresses all are greatest in the density minima be- 
tween the rings. The dense rings are 'quiet' locations where 
the turbulence is substantially weaker. Such an environment 
may be helpful for the growth of larger bodies. 
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- The rings within the active zone sometimes move about, 
leading to mergers. In contrast, the bump at the dead zone in- 
ner edge is stationary. Planetary embryos with masses in the 
type I migration regime can be retained a t the dead zone in- 
ner ed ge as proposed bv llda & Lirtd2008h ; ISchlaufman etal] 
d2009h . 

6. Outlook 

The causes of the magnetic field oscillations appearing in the 
butterfly diagram remain to be clarified. In particular, it is not 
known what processes drive the quasi-periodic reversals in the 
azimuthal magnetic fields shown in figure [4] 

While we have used a fixed magnetic diffusivity distribu- 
tion, the degree of ionization will in fact change as the disk 
evolves. Annuli of increasing surface density will absorb the 
ionizing stellar X-rays and interstellar cosmic rays further from 
the midplane, and will have higher recombination rates due to 
the greater densities. These effects will likely mean contrasts 
in the strength of the turbulence between the rings and inter- 
ring regions, even greater than those observed in our calcula- 
tions. Stronger magnetic fields and higher mass flow rates in 
the inter-ring regions could lead in turn to growth in the sur- 
face density contrast over time, analogous to the viscous insta- 
bility of the q-model in the radiation-pre ssure dominated regime 
dLightman & Eardlevll974tlPiranlll978l) . 

The resistivity can show another kind of time variation near 
the boundary of the thermally-ionized region. Changes in the 
temperature can alter the strength of the turbulence and thus the 
heating rate. In this way, radial heat transp ort can activate previ- 
ously dead gas CWi insch et al.1 '2005. 2006). Radial oscillations of 
the ionization front may be the consequence. Overall, due to the 
fixed magnetic diffusivity, it is probable that our models under- 
estimate the evolution resulting from the ionization thresholds. 

Owing to the effects of the radial boundary conditions, our 
global simulations have limitations for measuring quantities 
such as the accretion rate and the mean radial velocity. Fixing 
the rate of flow across the outer radial boundary is a promising 
avenue for further exploration in this direction. 
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